function [U,dUdTheta] = Payoff(AllStatesMap,MidPointsFixedStates,FixedStatesMap,SRDMidPoints,theta,y)
% Computes payoff U and partial of payoff wrt Theta for region r.

Nf  = size(FixedStatesMap,2);
Nd  = max(AllStatesMap(:,end));
Nsf = size(FixedStatesMap,1);

% State transition:
nstates = size(AllStatesMap,1);
ntheta  = length(theta);
V0 = (AllStatesMap(:,1) == 1);

% dUdTheta{i} stores the partial of U wrt Theta for payoff-action i.
% In our case, action 1 is inaction, 2 is plant or replant sugarcane and 3
% is abandon sugarcane.

U           = zeros(nstates,3);      % Initialize payoffs
dUdTheta{1} = zeros(nstates,ntheta); % Initialize payoffs der. w.r.t. theta.
dUdTheta{2} = zeros(nstates,ntheta);
dUdTheta{3} = zeros(nstates,ntheta);

gamma  = theta(1);
thetaR = theta(2);
psiE   = theta(3:7);
psiR   = theta(8);
psiA   = theta(9);
kappa  = theta(10);
beta   = theta(10+1:10+2*(Nf-3));

sprice = y(:,1); % Sugarprice vector
outret = y(:,2); % Outside option return vector

outretBig = repmat(outret,[Nsf*Nd,1]);
spriceBig = repmat(sprice,[Nsf*Nd,1]);

fixed_ind = FixedStatesMap(AllStatesMap(:,end-1),:);
fixed_cov = zeros(size(fixed_ind,1),size(fixed_ind,2)-1);
for cov = 1:size(fixed_ind,2)-1
    fixed_cov(:,cov) = MidPointsFixedStates(fixed_ind(:,cov),cov);
end

LUstate = fixed_ind(:,end);
SRDist    = SRDMidPoints(AllStatesMap(:,end),1);
SRDistFar = (AllStatesMap(:,end) == max(AllStatesMap(:,end)));

yield   = fixed_cov(:,2);
tc      = fixed_cov(:,1);
costvar = [fixed_cov(:,3:end) fixed_cov(:,3:end).^2];

U(V0,1) = thetaR*outretBig; 
sV0     = sum(V0);
dUdTheta{1}(V0,:) = [zeros(sV0,1) outretBig zeros(sV0,ntheta-2)];

U(V0,2) = - psiE(LUstate(V0)) - ~SRDistFar(V0).*SRDist(V0)*psiE(4) - SRDistFar(V0)*psiE(5); %PsiI

dUdTheta{2}(LUstate==1 & V0,3) = -1;
dUdTheta{2}(LUstate==2 & V0,4) = -1;
dUdTheta{2}(LUstate==3 & V0,5) = -1;
dUdTheta{2}(V0,[6 7]) = [- ~SRDistFar(V0).*SRDist(V0), - SRDistFar(V0)];

U(V0,3) = - Inf; % At other state, there is no outside option

maxage = double(max(AllStatesMap(:,1)));
for a=2:maxage
    Va  = (AllStatesMap(:,1) == a);
    yielda = yield(Va);
    tca    = tc(Va);
    costvara = costvar(Va,:);
    sVa = sum(Va);
    U(Va,1) = kappa*(gamma^(a-2))*(spriceBig-tca).*yielda + costvara*beta; % 
    dUdTheta{1}(Va,:) = [kappa*((a-2)*gamma^(a-3))*(spriceBig-tca).*yielda ...% partial wrt gamma
        zeros(sVa,8) ... % partial wrt thetaR and fixed costs ...
        (gamma^(a-2))*(spriceBig-tca).*yielda ...  % partial wrt kappa ...
        costvara]; % partial wrt beta
    
    U(Va,2) = - psiR;
    dUdTheta{2}(Va,8) = -1;
    
    U(Va,3) = - psiA;
    dUdTheta{3}(Va,9) = -1;
end


end